Lattice field theory simulations of graphene 
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We discuss the Monte Carlo method of simulating lattice field theories as a means of studying 
the low-energy effective theory of graphene. We also report on simulational results obtained using 
the Metropolis and Hybrid Monte Carlo methods for the chiral condensate, which is the order 
parameter for the semimetal-insulator transition in graphene, induced by the Coulomb interaction 
between the massless electronic quasiparticles. The critical coupling and the associated exponents 
of this transition are determined by means of the logarithmic derivative of the chiral condensate 
and an equation-of-state analysis. A thorough discussion of finite-size effects is given, along with 
several tests of our calculational framework. These results strengthen the case for an insulating 
phase in suspended graphene, and indicate that the semimetal-insulator transition is likely to be of 
second order, though exhibiting neither classical critical exponents, nor the predicted phenomenon 
of Miransky scaling. 

PACS numbers: 73.63.Bd, 71.30.+h, 05.10.Ln 
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I. INTRODUCTION 

The recent experimental isolation of single atomic lay- 
ers of graphite, known as graphene, has provided physi- 
cists with a novel opportunity to study a strongly cou- 
pled system with remarkable many-body and electronic 
properties, which at the same time can be easily ma- 
nipulated experimentally P, Even more recently, 
the advent of experiments utilizing samples of suspended 
graphene, free from the interference of an underlying sub- 
strate 0, has provided unprecedented insight into the 
intrinsic properties of graphene. Among other remark- 
able discoveries, suspended graphene has been shown to 
possess a very high carrier mobility even at room tem- 
perature, as well as a markedly non-metallic behavior of 
the conductivity at low temperatures. 

A central property of graphene is that the low-energy 
electronic spectrum can be described in terms of two fla- 
vors of massless, four-component fermionic quasiparticles 
with linear dispersion [4(. Indeed, due to the hexago- 
nal honeycomb arrangement of the carbon atoms in the 
graphene lattice, the band structure of graphene exhibits 
two inequivalent (but degenerate) "Dirac cones" where 
the conduction and valence bands touch, as illustrated in 
Fig. 
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Since the energy-momentum relation around 
a Dirac point is linear as in relativistic theories, the low- 
energy description of graphene bears a certain resem- 
blance to massless Quantum Electrodynamics (QED). 
Nevertheless, an important difference is that the Fermi 
velocity of the quasiparticles in graphene is as low as 
v ~ c/300, whereby the electromagnetic interaction is 
rendered essentially instantaneous. 

Such a description is well-known to account for the 
physics of graphene on a substrate, where the system ex- 
hibits semimetallic properties due to the absence of a gap 
in the electronic spectrum. While suspended graphene 
has recently come under intense experimental investiga- 
tion [U , its spectrum is yet to be computed in a controlled 



fashion. From the theoretical perspective, the challeng- 
ing feature of suspended graphene lies in the fact that 
the Coulomb interaction between the quasiparticles is un- 
screened which, in conjunction with the small Fermi ve- 
locity, results in a graphene analogue of the fine-structure 
constant a g > 1. At such strong coupling, a dynami- 
cal transition into a phase fundamentally different from 
the weakly-coupled semimetallic phase of graphene is a 
strong possibility. In graphene sheets deposited on a sub- 
strate, such a transition is effectively inhibited due to the 
screening of the Coulomb interaction by the dielectric. 

Our recent work in Ref. 5j has demonstrated that 
graphene is expected to undergo a semimetal-insulator 
transition when the substrate is removed. More specif- 
ically, evidence was found that the low-energy effective 
theory of graphene undergoes a phase transition involv- 
ing spontaneous chiral symmetry breaking, which takes 





FIG. 1: (Color online) a) Dirac cone, joining the upper (red) 
conduction band and the lower (blue) valence band, b) The 
hexagonal arrangement of carbon atoms in graphene, with 
sublattices A (red dots, thin dashed lines) and B (green dots, 
thin solid lines). 
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place at a critical coupling of (3 C = 0.072 ± 0.005, and 
within the accuracy of that work the transition appeared 
to be consistent with classical mean-field exponents. The 
results reported in Ref. @ are based on the numerical 
Monte Carlo simulation of a discretized lattice formula- 
tion of the low-energy effective theory of graphene, and 
the calculation of the chiral condensate, which is the or- 
der parameter for excitonic gap formation. In one possi- 
ble realization of such an insulating state, the equivalence 
of the triangular sublattices A and B, shown in Fig. l(b)| 
is broken by the accumulation of charge carriers of oppo- 
site sign on the respective sublattices. 

Our results should be compared with those of Refs. [6j, 
0] which are based on a gap equation, where a semimetal- 
insulator transition was found at critical couplings of 
j3 c ~ 0.06 and f3 c ~ 0.03 respectively. While the re- 
sult of Ref. @ is within the physical range of Coulomb 
couplings, that of Ref. 0] is slightly above the largest 
conceivable value of a g ~ 2.16, which corresponds to 
graphene in vacuum. On the other hand, Refs. [j| 
employed an expansion in the inverse number of fermion 
flavors Nf, and found that at large 7V^, the Coulomb in- 
teraction between the quasiparticles becomes irrelevant 
and therefore unable to induce a gap in the electronic 
spectrum. 

In this paper, we explain the details of our Lattice 
Monte Carlo method, which to our knowledge has not 
been applied to the low-energy theory of graphene (how- 
ever, Ref. has considered a theory related to the 
strong-coupling limit). We also present new calculations 
supporting the conclusions of Ref. @ , but extending the 
previous data set to much larger lattices. In Section HT1 
we discuss the low-energy effective theory of graphene, 
the corresponding partition function, and the compu- 
tation of observables upon integration of the fermionic 
degrees of freedom. In Section IIII1 we describe the dis- 
cretization of the effective theory and introduce a lattice 
formulation that respects gauge invariance and avoids 
the fermion doubling problem while maintaining a cer- 
tain degree of chiral symmetry at finite lattice spacing. 
In SectionfV] the results of our simulations are presented, 
with emphasis on the chiral condensate and susceptibil- 
ity, including a determination of the critical coupling for 
the semimetal-insulator phase transition, and the con- 
sequences of our results for the corresponding critical 
exponents. In Section IVI1 we outline the various tests 
and cross-checks we have performed in order to validate 
our results. In Section IVIII we discuss the possibility of 
observing the transition experimentally. Finally, in Sec- 
tion I Villi we summarize our findings and present a case 
for continued study. 



II. LOW-ENERGY EFFECTIVE THEORY 

The electronic band structure of graphene close to the 
Fermi level forms the basis of the low-energy effective 
theory of graphene. This band structure is a reflection of 



the hexagonal arrangement of the carbon atoms as shown 
in Fig. 1 (b) and can be well described by a tight-binding 



model of the form 

H = 1 E (4a, +h. 

E (<i«a,i+^^+H.c), (2.1) 

«*,3»,o-=t.J, 

as first done by Wallace in Ref. The operators 

a ti( a (7.i) an d f>ai(b ai ) create (annihilate) an electron of 
spin a at location i on the A and B sublattices, respec- 
tively [see Fig. l(b)| . The first term (involving t) takes 
into account nearest-neighbor interactions, and the sec- 
ond term (involving t') the next-to- nearest neighbor ones. 
Both terms account for all spin states. The hopping pa- 
rameters that give an optimal fit to the experimentally 
determined band structure of graphene are t ~ 2.8 eV 
and t' ~ 0.1 eV [ljj]. Third-nearest neighbors have also 
been considered in Ref. [T^, yielding an additional hop- 
ping amplitude of t" ~ 0.07 eV. 

We shall follow a somewhat different route based on an 
Effective Field Theory (EFT) treatment of graphene !). 
[l]|, which has the advantage of describing the physics 
of graphene directly in terms of the relevant low-energy 
degrees of freedom, namely charged massless fermionic 
quasiparticles. The EFT description of graphene has an 
additional advantage as it allows for the direct study of 
effects due to the unscreened, long-range Coulomb in- 
teractions between the quasiparticles. In what follows, 
we shall formulate a continuum Lagrangian field theory 
that should be thought of as valid only at low momenta, 
much smaller than the inverse of the interatomic distance 
in graphene, which is ~ 1.42 A. 



A. Continuum formulation 

In the EFT framework, graphene is described by a the- 
ory of Nj: Dirac flavors interacting via an instantaneous 
Coulomb interaction. The action (in Euclidean space- 
time) of this theory is 



f 

S E = d 2 xdti, a D[A ]ij a 

d 3 xdt(diA ) 2 , (2.2) 



a=l 
I 

+ V 

where = 2 for graphene monolayers, g 2 — e 2 /e 
for graphene in vacuum (suspended graphene), ip a is a 
four-component Dirac field in 2+1 dimensions, A is a 
Coulomb field in 3+1 dimensions, and 

D[A ] = 7o(0o+*4>)+«7i0i, * = M (2.3) 

where the Dirac matrices 7„ satisfy the Euclidean Clif- 
ford algebra {7^, 7^} = 25^. The four-component spinor 
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structure accounts for quasiparticle excitations of sub- 
lattices A and B around the two Dirac points in the 
band structure 0, 0| ■ The two Dirac points are identi- 
fied with the two inequivalent representations (with op- 
posite parity) of the Dirac matrices in 2+1 dimensions. 
In graphene monolayers, Nf = 2 owing to electronic spin, 
while Nf = 4 is related to the case of two decoupled 
graphene layers, interacting solely via the Coulomb in- 
teraction. Consideration of arbitrary Nf is also useful, 
given that an analytic treatment [sj is possible in the 
limit Nf — > oo. 

The strength of the Coulomb interaction is controlled 
by a g — e 2 /(47rue ), which is the graphene analogue 
of the fine-structure constant a ~ 1/137 of QED. It is 
straightforward to show that a g is the only parameter, 
by rescaling according to 



t' = vt, 
A' = AJv. 



(2.4) 



The action (|2.2p is invariant under spatially uniform 
gauge transformations (see Sec. IIII A|) . Notice that since 
the A field is 3+1 dimensional, one recovers the four- 
fermion Coulomb interaction 



operate in flavor space. In this way, one obtains a set 
of precisely sixteen linearly independent Hermitian ma- 
trices, forming the generators of U(4). Significantly, this 
chiral symmetry can be spontaneously broken down to 
U(2)xU(2), in which case the excitonic condensate (4'ip) 
acquires a non- vanishing value, signaling the formation of 
quasiparticle- hole bound states. The same group struc- 
ture is obtained by adding to Eq. (|2.2p a parity invariant 
(Dirac) mass term 



d 2 xdtm il> a ip a , 



(2.8) 



(2.5) 



which breaks chiral symmetry explicitly. The remaining 
unbroken generators are then {l,cr 3 }, which correspond 
to uniform phase rotations of both flavors with the same 
phase, and with equal and opposite phases, respectively. 
For the extended theory with Nf flavors, the symmetry 
breaking pattern is V(2N f ) -> \J (N f ) x\J (N f ) . 

Other symmetry breaking patterns, particularly in- 
volving the possibility of magnetic as well as Cooper-like 
pairing instabilities, have been investigated in Refs. [jj 



B. Effective action and probability measure 



by integrating out A Q . Nevertheless, for our purposes the 
original form of the action (quadratic in the fermions) as 
given in Eq. (|2.2[) is preferable. 

A central property of the low-energy EFT is that 
Eq. (|2.2[) respects a global U(2iV^) chiral symmetry un- 
der the transformations 



exp(«r -a ) tp a 



(2.6) 



where the matrices 1^ are the (2Nf) 2 Hermitian gen- 
erators of \J(2Nf), such that for the case of graphene 
monolayers, the group is U(4). The generators can be 
constructed by first choosing a representation for the 7^, 
such as 



7o 



a 3 
-cr. 



7i 



a l 

-cr- 



(2.7) 



where the a i are Pauli matrices. Adding the identity to 
this set yields the generators of U(2), since they form a 
set of four linearly independent Hermitian matrices. It 
should be noted that the choice of any particular repre- 
sentation for the 7^ is completely arbitrary and is not 
necessary for any calculational purpose, as all relevant 
information is provided by the Clifford algebra. How- 
ever, the identification of the spinor degrees of freedom 
with any particular Dirac point and graphene sublattice 
is dependent on the chosen representation. 

In order to arrive at the generators of U(4), one can 
take the direct product of each of the abovementioned 
generators of U(2) by {1, cr 1 , cr 2 , cr 3 }, where the latter 



The partition function corresponding to Eq. (12. 2|) is 
given by 

Z = ( VA Q V^Vi> eMS E $a^a,A a }), (2.9) 



where it is possible to integrate out the fermionic degrees 
of freedom, as S E is quadratic in the ip a . We thus obtain 

Z = JvA cxp(-S|[A ]) dct(D[A }) N f , (2.10) 



where 



(2.11) 



is the pure gauge part of the action. It is of central impor- 
tance for the convergence of the Monte Carlo algorithm 
that the above determinant has a definite sign, indepen- 
dently of any particular configuration of the gauge field 
A . One way to prove that this property is satisfied is 
to choose a specific representation of the Dirac matrices, 
such as Eq. (|2.7|) . in terms of which -D[A ] can be written 
as 



D[A ] 



M[A ) 
-M[A ] 



M[A ] 
Aft [Aq] 



where 



(2.12) 



M[A ]=<T (d +iA )+va i d i , i = 1,2, (2.13) 
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and use the facts that A is real, and that the Pauli ma- 
trices and the momentum operator are Hermitian. The 
latter implies = — 9 , and therefore 

det(D) = det(M)det(M t ) = |det(M)| 2 > (2.14) 

which, furthermore, is not affected by the introduction of 
a parity invariant mass term such as Eq. (|2.8|) . However, 
the positivity of det(-D) breaks down in the presence of a 
chemical potential, which can be thought of as a uniform, 
imaginary contribution to the A field. 

The fact that det(D) is positive definite allows for the 
definition of an effective gauge action that is purely real, 
given by 

S cS [A Q ] = -N f lndet(D[A }) + S g E [A ], (2.15) 
so that the partition function becomes 



Z = JvA Q cxp(-5 off [A ]), 



(2.16) 



where -P[^4 ] = exp(— S cS [A ]) > can be interpreted as 
a positive definite probability measure for a Monte Carlo 
calculation, as outlined in Section HTT1 



C. Operator expectation values 

The expectation value of a given operator Ofy, ip] de- 
pendent on the fermion fields can be calculated by taking 
functional derivatives of the generating functional 



Z[fj,ij\ = J VA V4>V^ 

x expt-S^LV^V^,??]), (2-17) 

where source terms have been added to the original action 
according to 

S E [ A o^^^,v] = S E [A ,$,ip]+ / d 2 xdt{ipr] + H.c), 



(2.18) 

such that the modified effective gauge action is a func- 
tional of A as well as of the sources rj, fj. It is again 
possible to integrate out the fermionic degrees of free- 
dom and take functional derivatives with respect to the 
sources in the resulting expression, 

Z[fj,rj\ ex J VA exp(-S eS [A ]) 

x exp(- / d 2 xdtfjD- 1 [A ]ri ] , (2.19) 



which makes it possible to obtain expectation values in 
terms of a path integral over A only. While this pro- 
cedure is completely general, it is possible to employ a 
slightly different approach in order to facilitate the com- 
putation of the chiral condensate and susceptibility. 



The chiral condensate a, which is the order parameter 
of the semimetal-insulator phase transition in graphene, 
is defined by 



bh 



(2.20) 



where the fermion fields are evaluated at the same space- 
time point. It is useful to note that the mass m plays 
the role of a source, coupled to The expectation 

value of this operator can therefore be obtained by first 
differentiating the partition function with respect to m Q 
and dividing by the volume, giving 

a = JL [ VA^V^ I dx$ b (x)i> b (x) exp(-S B ) 



1 dlnZ 
V dm 



(2.21) 



where we have used the fact that space is homogeneous 
and therefore the volume average of l tp b (x)tjj b (x) can be 
replaced by its value at an arbitrary point x. On the 
other hand, once the fermions have been integrated out, 
the derivative with respect to m yields 



X - J VA.TiiD-^A,]) eMS cS [A ]) 



vz 

^(TrCD-Vo])) 



where the identities 



det(0[A]) = exp(Tr(log(£[A])), 



ddet(D[X\) 



det(D[X])Ti[D- 1 [\] 



3D 



(2.22) 

(2.23) 
(2.24) 



have been used. The chiral susceptibility \i may be found 
by taking one more derivative with respect to m , giving 



Xi 



(2.25) 



V 



(Tr 2 ^- 1 )) - (Tr(D- 2 )> - (Tr^- 1 )) 1 



which is expected to diverge at a second-order phase tran- 
sition, and may also yield constraining information on the 
universal critical exponents of the transition. 



III. GRAPHENE ON THE LATTICE 

In this section we formulate the lattice version of 
Eq. (|2.2p . We begin by discretizing the pure gauge sector, 
where the requirement of gauge invariance implies the 
use of "link variables" to represent the gauge degrees of 
freedom. The "staggered" discretization of the fermionic 
sector is then outlined, as it is the preferred choice to 
represent fermions with chiral symmetry at finite lattice 
spacing. 
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A. Gauge invariance and link variables 

Recall that the pure gauge part of the Euclidean action 
is given by 



C<9 



1 

V 



(fxdt{d i A ) 



(3.1) 



which can be thought of as the non-relativistic limit of the 



Lorentz-invariant form \F^ LV F^ V where F^ u 
d v A , such that 



F F lJ 



F l0 F 



= 2F ,F°' = 2(5,A ) 2 



(3.2) 



where we have used F^ = (no magnetic field) and 
d Aj — (no electric field induction by a magnetic field), 
valid in the non-relativistic limit (v -C c). Thus, for 
graphene the only non- vanishing contribution is the elec- 



tric field E 4 



-d jA , which represents the instanta- 



neous Coulomb interaction between the quasiparticles. 

The action (|3.ip is invariant under the time-dependent, 
spatially uniform gauge transformations 



In terms of the gauge link variables and plaquettes, 
the discretized gauge action corresponding to Eq. (|3.ip 
is given by [151 ] 



l 

l -2 



U 



(3.8) 



where — l/g 2 , such that (3 — > v/g 2 when the rescal- 
ing of Eq. ([21} is applied. In Eq. (O, the only 
non-vanishing contributions arise from the terms with 
( M , v) = (1, 0); (2, 0); (3, 0); (2, 1); (3, 1) and (3, 2). Equa- 
tion (|3.8p may be simplified to 



S 9 e,c 



^E 



3 

E 

i=l 



n v n+e 



(3.9) 



where 3?(x) denotes the real part of x. Equation (|3.9p is 
referred to as the compact formulation of the discretized 
gauge action. This formulation is known [161 ] to be sub- 
optimal, as compared to the non-compact formulation, 
for lattice simulations of QED and related theories. How- 
ever, the non-compact formulation may be obtained from 
Eq. (|3.9p by expanding ^(U n U^ +e ) to second order in 9, 



A -> A + a(t), 
ip — ► exp j dt'a(t') J> / '. 



(3.3) 



where a(t) is a function of time only. Thus, in spite of 
its apparent simplicity, the effective theory of graphene 
possesses a truly local gauge invariance, which should be 
respected by the lattice action. To this end, one intro- 
duces temporal link variables 



£Vn = U a = exp (i0 n ) 



(3.4) 



where 8 n is the dimensionless lattice gauge field evaluated 
at the lattice point n = (no, rii, n^, 713). The spatial link 
variables 



Ui, n = 1, 



(3.5) 



are set to unity. It is convenient to express the discretized 
version of Eq. (|3.ip in terms of "plaquette" variables, 
defined by 



u, 



(3.6) 



where, in the present case of a pure Coulomb interaction, 
the only non-trivial components are U oi and U i0 . Those 
plaquette components then correspond to the discretized 
formulation of the electric field. The remaining compo- 
nents corresponding to the magnetic field are equal to 
unity. These statements can be summarized in the ex- 
pression 



+ S tl0 S 1/0 + S tli S I/j . 



(3.7) 



K(^+e..) = l-2 V"" -e 



(3.10) 



whereupon the non-compact lattice gauge action is given 
by 



08 _ P ST 

e,n - 2 2-^2.^ 



n+e ? - n 



(3.11) 



Here, and throughout the rest of this paper, we have set 
the lattice spacing to equal unity, and it is thus dropped 
from all expressions. All dimensionful quantities should 
therefore be regarded as expressed in units of the lattice 
spacing. 



B. Staggered fermions 

While the discretization of the gauge sector is relatively 
straightforward, the inclusion of dynamical fermions on 
the lattice is a notoriously difficult problem. One of 
the main issues when simulating fermions on the lattice 
is the so-called doubling problem (for an overview, see 
Ref. [HI], Chapter 4). This problem is related to the chi- 
ral invariance of the fermionic sector, and arises due to 
the appearance of multiple (unwanted) zeros in the in- 
verse propagator. In other words, one is simulating more 
fermion flavors than expected, the exact number being 
dependent on the dimensionality of the theory. There 
exists a number of ways to avoid the doubling problem, 
but all of them break chiral invariance in one way or the 
other, an inevitable fact encoded in the Nielsen-Ninomiya 
theorem [l7| . The solution we have chosen for our simu- 
lations of graphene is the "staggered" fermion represen- 
tation of Ref. [l8| . This choice is optimal for the study 
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of spontaneous chiral symmetry breaking in graphene, as 
it yields the correct number of degrees of freedom while 
also partially preserving the original chiral symmetry of 
the theory, as will be shown in this section. 

In order to discretize the fermionic sector of Eq. (|2.2[) 
in a way amenable to computer simulations, there are 
a number of choices that need to be made. As a first 
step, the fermions are integrated out, and the problem is 
formulated using the partition function written purely in 
terms of the gauge field [Eq. (|2.16[) ]. The fermions are 
then represented exclusively through the determinant of 
the Dirac operator D. One can attempt to compute the 
determinant exactly for a given 9 configuration, which 
is feasible due to the low dimensionality of the problem, 
and is what we have done for part of our calculations. 
Alternatively one may use the so-called pseudofermion 
method, which we will briefly explain in the next section. 

In order to arrive at the staggered fermion formulation, 
a useful starting point is the "naively" discretized action 

S f E [^,e] = -J2AnD n , m [9]\ m , (3.12) 

n.m 

where 

D n, m [ 6 \ = ^70(<5 n+eoim t/ n _ S n-e ,m U L) 

i 

+ m (S n,m> ( 3 - 13 ) 

with U n = exp(i8 n ). It should be noted that for small 
m , Eq. (I3.13P becomes ill-conditioned, such that the 
"chiral limit" m a — ► has to be reached by extrapo- 
lation. The boundary conditions of the fermion fields are 
periodic in the spatial directions and anti-periodic in the 
temporal direction. It is possible, using a local unitary 
transformation on the fermion fields, to simultaneously 
diagonalize the Dirac matrices in Eq. (|3.13[) and thereby 
decouple the spinor components. This procedure, known 
as the Kawamoto-Smit transformation (l9l | or simply as 
"spin-diagonalization" , is defined by 

Ai -> T nXn, 

^„ - X a Tl (3.14) 

which in the Dirac operator (|3.13|) effects the transfor- 
mation 

y-Ttyr^ (3 . 15) 

on the Dirac matrices 7**. The transformed fermion fields 
X n are referred to as staggered spinors. It is straightfor- 
ward to show that the choice T n — Jq ^ 1 ^ 2 satisfies 

TlrT n+e ^ <1, (3.16) 
where the Kawamoto-Smit phases are given by 

€ = 1, 

vl = {-l) n « +n K (3.17) 



In this fashion the Dirac structure is removed, result- 
ing in a sum of four identical terms in the action, one 
for each component of the original four-component Dirac 
spinor ip a . These copies are referred to as staggered fla- 
vors. It has been shown in Ref. 20] that for each stag- 
gered flavor one recovers, in the continuum limit, two 
four-component Dirac flavors. Thus, by retaining one 
staggered flavor, it is possible to have exactly eight con- 
tinuum fermionic degrees of freedom, which is the correct 
number for graphene. The action of a single staggered 
flavor is given by 

s f E [x,x,e] = -J2x n K n ^[e} Xm , (3.18) 

n,m 

where the staggered Dirac operator is 

^n,m[^] = 2 (^ n + e ' m ~~ ^n-e ,m ^im) 

+ 2 / y (^n+e; ,m — ^n-e^m) 

i 

+ ™0 (5 n,m- (3-19) 

The operator K thus replaces D in all expressions for 
the probability, chiral condensate and susceptibility that 
were derived in the previous sections. As expected from 
the Nielsen-Ninomiya theorem, the staggered lattice ac- 
tion does not retain the full U(4) chiral symmetry of 
the original graphene action at finite lattice spacing. As 
shown in Ref. [2(|, only a subgroup U(l)xU(l) remains 
upon discretization. Spontaneous condensation of xXi or 
equivalently the introduction of a parity invariant mass 
term, reduces this symmetry to U(l). The focus of this 
work is on the phase transition associated with such a 
chiral symmetry-breaking pattern. 

Finally, it should be pointed out that the situation 
concerning graphene is unusually favorable, in the sense 
that the staggered formalism somewhat fortuitously pro- 
vides the correct number of fermionic degrees of freedom 
as N f — 2 for graphene monolayers. In general, stag- 
gered fermions provide only a compromise solution in the 
sense that some degree of chiral symmetry is preserved, 
at the price of retaining some of the doubling originally 
present in the discretized fermion action. Indeed, if the 
case of Nj — 1 were to be simulated, it would be neces- 
sary to resort to the uncontrolled and controversial "root- 
ing" trick [2ll ] , whereby the desired number of continuum 
flavors is restored by taking the appropriate root of the 
Dirac operator. 



C. Computation of observables 

The computation of a and \i from ensembles of gauge 
field configurations necessitates, in principle, the full in- 
version of K and K 2 . Such a procedure may potentially 
become extremely time-consuming for large lattices. In 
this respect, a choice exists between direct sparse solvers, 
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such as PARDISO [22j], and iterative solvers of the con- 
jugate gradient type, such as BiCGStab (23[. As the 
lattice size is increased, the performance of the direct 
solver scales much worse than the iterative solver, by a 
factor roughly proportional to the lattice volume. Nev- 
ertheless, the direct sparse solvers remain an attractive 
choice for a number of reasons: the performance of a di- 
rect solver is independent of the condition number of K, 
which is the ratio of its largest and smallest eigenvalues, 
and this is particularly significant close to a transition 
and for small m . Furthermore, direct solvers feature 
optimized parallelization and are efficient at handling in- 
version problems with multiple right-hand-sides. In view 
of this, PARDISO has been found to be the solver of 
choice for the efficient computation of observables on the 
presently used lattice volumes. 

Regardless of the type of solver used, the full inversion 
of K quickly becomes impractically expensive when the 
lattice size is increased. In this situation, it is possible 
to resort to a stochastic estimator [24[ , which constitutes 
an alternative to the exact calculation of Tr(X _1 ). A 
suitable stochastic estimator for a is given by 

II III 

where the £ n are random Gaussian variables which sat- 
isfy ((£ n )) = an d ((£mO) = ^m.n; where the double 
bracket notation indicates an average over £ n . 

For a given gauge configuration, averaging Eq. (|3.20|) 
over £ n yields Tr(if _1 ), which only requires application 
of the inverse to a limited number of random Gaussian 
vectors. With this approach it is also straightforward to 
compute Tr(X -2 ), by simply applying the inverse to each 
random vector one more time. Adequate accuracy for a 
and Xi is achieved using ~ 100 random vectors for each 
gauge configuration, independently of the lattice volume 
used. 



IV. MONTE CARLO STRATEGIES 

This section presents the two Monte Carlo algorithms 
that we have used to study the discretized low-energy 
effective theory of graphene. We begin by outlining the 
Metropolis Monte Carlo algorithm which, although con- 
ceptually simpler, becomes computationally inefficient 
beyond a certain lattice volume, after which we proceed 
to describe the more advanced and highly efficient ap- 
proach involving the Hybrid Monte Carlo (HMC) algo- 
rithm with pseudofcrmions. 

A. Metropolis Monte Carlo 

As shown in Section III Bl the structure of the fermion 
determinant allows for a positive definite probability 
measure. Indeed, as shown in Section III C[ an effective 



action can be defined such that expectation values of ob- 
servables can be written as averages over field configura- 
tions weighted by 

P[9] ee ex P (-S cff [0]) = det(K[e}) exp(-Sf [0]), (4.1) 

where the matrix K corresponds to the staggered Dirac 
operator of Eq. (|3.19|) . In the Metropolis algorithm [25| . 
a given gauge field configuration 9 is updated by the in- 
troduction of a small change at a randomly chosen lattice 
site. The updated configuration 9' is then accepted with 
probability 

P\9'] 

V ee -J^=cxp(-AS), 

AS = S eS [9']-S eS [9\. (4.2) 

If the new configuration 9' is rejected, 9 is retained, and a 
new change proposed. In this fashion, a so-called Markov 
chain of gauge configurations is generated, in which the 
samples are distributed according to the desired probabil- 
ity measure. After an appropriate number of thermaliza- 
tion steps, gauge configurations can be saved at regular 
intervals, which should allow for adequate decorrelation. 
The central limit theorem then guarantees that for Af un- 
corrected samples, the statistical uncertainties will de- 
crease as 1/VJ7. The decorrelation can be measured in 
terms of the number of full sweeps of the lattice required 
between two consecutive observations, in order for the 
autocorrelation of the ensemble of gauge configurations 
to become insignificant. For the Metropolis algorithm, a 
proper balance between update size and decorrelation is 
achieved for acceptance rates of ~ 60 — 70%. 

In spite of its simplicity, the Metropolis approach has 
several inherent disadvantages. The most serious one 
arises as the fermion action is non-local, in the sense 
that updating a single lattice site requires a full recal- 
culation of det(K). This disadvantage is exacerbated by 
the fact that decorrelation is dependent on the number 
of full sweeps of the lattice, and the number of sites to be 
updated increases as the lattice size is increased. Even 
with highly efficient parallel sparse solvers, the execution 
time scales as ~ V 3 , such that it is bound to become 
impractical above a certain maximum lattice size. Also, 
as the updates in the Metropolis algorithm are entirely 
random, it is usually only possible to update very few lat- 
tice sites at once without ruining the acceptance rate. In 
the Sec lIV Bl we give an overview of the HMC algorithm, 
which is designed to overcome these difficulties. 

B. Hybrid Monte Carlo 

The problem of efficient updating of the gauge field 
in theories with dynamical fermions has been addressed 
in Ref. [26| where the Hybrid Monte Carlo (HMC) al- 
gorithm was introduced. In this approach, the gauge 
field is evolved deterministically along a Molecular Dy- 
namics (MD) trajectory, such that the entire lattice is 
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updated at once. Thus, the number of updates required 
for decorrelation within the HMC algorithm is dramat- 
ically reduced, although the number of MD trajectories 
required for decorrelation roughly equals the number of 
sweeps necessary in the Metropolis approach. 

The basic idea of the HMC algorithm is to evolve a 
given initial configuration 8 n in a fictitious time r accord- 
ing to the classical equations of motion, with a Hamilto- 
nian given by 



H 



S E [9] 



(4.3) 



where iSs[0] is the Euclidean action to be sampled, and 
7r n is a momentum conjugate to 9 n . This momentum is 
introduced as an auxiliary field, with the sole purpose 
of defining the above dynamics. The field 7r n is of no 
consequence to the path integral that defines the theory, 
as its contribution factors out completely. It has been 
shown in Ref. [26j that the procedure of classically evolv- 
ing (8 n , 7r n ) — » (6^, using the above Hamiltonian, and 
choosing the initial 7r n from a random Gaussian distri- 
bution, produces a Markov chain of gauge field configu- 
rations distributed according to the desired probability 
measure. 



Because the MD evolution is in principle exact, a tra- 
jectory that is long enough should provide the desired 
decorrelation between consecutive samples, provided that 
the pseudofermion field is refreshed at regular intervals. 
Ideally, a 100% acceptance rate should thus be achiev- 
able. In practice, however, the MD evolution is imple- 
mented with a finite time step At, which introduces a 
systematic error. However, as long as the evolution re- 
mains reversible, the effects of that error on the distri- 
bution of gauge field configurations can be eliminated 
by means of a Metropolis step, comparing the initial 
and final configurations after each MD evolution, where 
Eq. (|4.3p plays the role of the effective action in Eq. (|4.2[) . 



While the HMC algorithm achieves very efficient up- 
dating of the gauge field, a potentially serious draw- 
back is that the updating procedure requires (in prin- 
ciple) the full evaluation of K^ 1 which is computation- 
ally prohibitively expensive, even more so than det(K). 
Because of this, a number of methods have been devel- 
oped that seek to circumvent the necessity of calculating 
A' -1 . In one of these, the so-called R-algorithm 
the inverse is approximated by a stochastical estimator 
which, however, introduces a systematical error due to 
the loss of reversibility. Arguably, the method of choice 
is the ^-algorithm [2jJ, which reduces the MD evolution 
into a sparse operation by re-expressing the square of 
the fermion determinant as a path integral over complex 
scalar fields known as pseudofermions, while simultane- 
ously maintaining the desirable features of the HMC ap- 
proach. 



C. Pseudofermions 

As the pseudofermion method is explained in great de- 
tail elsewhere (for pedagogical reviews, see Refs. dEHH) 
we shall only concern ourselves with outlining the basic 
idea, which is based on the identity 



det(Q) oc / V$V<\> exp( S^), 



(4.4) 



where the constant of proportionality is of no conse- 
quence. Here, (j>, (ft are pseudofermion fields (which are 
bosonic but nevertheless satisfy anti-periodic boundary 
conditions in the temporal direction), Q = K and the 
pseudofermion action is 



(4.5) 



where £ follows a Gaussian distribution, related to the 
pseudofermion field by <f> — K'£. 

In order to simulate graphene, one requires det(K), not 
det(Q) = det(K'K). Thus, using the pseudofermions ac- 
cording to the above prescription effectively doubles the 
number of degrees of freedom. Fortunately, the staggered 
fermion action allows for an odd-even decomposition [28| , 
such that a single staggered flavor can be simulated. In 
the odd-even decomposition, the lattice is separated into 
sublattices of even and odd sites, according to the sign of 
(_l)»o+ n i+ n 2. Thus, as the derivative operator connects 
odd (even) sites with even (odd) ones, while the mass 
term connects odd (even) sites with odd (even) ones, the 
following odd-even decomposed form results: 



_ [ m o ^oe 
K eo m 



and therefore 



eo^oe ~ "'0 







kLk 



oe Mv eo "'0 



(4.6) 



(4.7) 



which, using the fact that K\ e = —K eo , has been factor- 
ized into blocks of even-even and odd-odd elements. As 



a consequence, 



det(Q) = det [K\ K 



(4.8) 



Thus, in order to recover det (if), it suffices to retain only 
the even-even (or odd-odd) block of Q. In practice, this 
is implemented simply by discarding either the odd (or 
even) elements of (f>. 

In the presence of pseudofermions, the MD Hamilto- 
nian becomes 

(4.9) 



and the equations of motion are 
• SH 

SH 

= 'so; 



F 9 + F p 

x n ~ n ) 



(4-10) 
(4.11) 
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where the "force term" associated with the gauge action 
takes the form 



S_S% 
S9„ 



(4.12) 



,9 Z 



6 n — ^n+ej + ^n-e^. 



where is the imaginary part of x; the second line 

in this equation corresponds to the compact formulation 
and the last line, obtained by expanding in powers of 6, 
shows the result for the non-compact case. The pseud- 
ofermion contribution is given by 



F p = - 



SSI 
S9„ 



(4.13) 



E 



S9„ 



The essence of the $- algorithm is the treatment of <fi as 
a constant background field throughout each MD trajec- 
tory. After each MD evolution, the pseudofermion field 
is refreshed using random Gaussian noise according to 
<f> = K'£. Computationally, the great advantage of this 
algorithm is that in each step At, the calculational ef- 
fort is reduced to applying the inverse of K^K to a single 
vector (f>, which is significantly less expensive than com- 
puting the full inverse. 

The numerical integration of the MD equations of mo- 
tion requires a reversible method, and the usual choice is 
the leap-frog integration formula [2(| which is also area- 
preserving. The calculation of the pseudofermion force 
in Eq. (|4.13[) is preferentially accomplished using an iter- 
ative solver such as BiCGStab (23[, in which case the al- 
gorithm scales roughly as ~ V. Nevertheless, in practical 
calculations the scaling is inevitably somewhat worse, as 
the truncation error of the leap-frog method tends to in- 
crease with increasing lattice size, necessitating a smaller 
timestep At. 

In the present study of the low-energy effective theory 
of graphene, we have used both the Metropolis and HMC 
algorithms, verifying that for any given set of parameters 
the results agree within statistical uncertainties. We now 
turn to a presentation of our simulation results. 



V. RESULTS 

In our simulations, the fermions live in a volume of 
extent V = L\ x L t , while the gauge bosons also prop- 
agate in the z-direction of length L z . Increasing L z be- 
yond 8 was found to have no discernible effects. The re- 
sults will thus be referred to by the short-hand notation 



L 2 x L t . Also, the action (|2.2[) has been rescaled accord- 
ing to Eq. (|2.4p . such that (3 = v/g 2 and v = 1 in the 
staggered Dirac operator of Eq. (|3.19[) . Our simulations 
have been performed at finite (but small) values of m , 
such that the limit m — * is reached by extrapolation. 

We have performed simulations on lattice sizes up to 
20 2 x 20 using the Metropolis method and 28 2 x 28 us- 
ing HMC. The former method scales roughly as V 3 and 
therefore quickly becomes uneconomical when the lat- 
tice volume is increased. However, an advantage of the 
Metropolis method is that the speed of the algorithm is 
independent of the condition number of the staggered 
Dirac operator K, as the fermionic determinant is eval- 
uated using a direct solver. In contrast, the HMC algo- 
rithm with pseudofermions scales roughly as ~ V, if used 
together with an iterative solver such as BiCGStab [23| . 
However, the HMC algorithm then becomes sensitive to 
the condition number of K, such that obtaining data be- 
comes more difficult at small bare fermion masses or close 
to the critical coupling. This problem can be somewhat 
alleviated using a direct solver such as PARDISO [22] |. 
but in that case the HMC algorithm scales roughly as 

Within the Metropolis approach, ~ 240 uncorrelated 
configurations were generated for each value of (/?, m ). 
When using the HMC algorithm, a similar number of 
MD trajectories were generated for each datapoint. The 
optimal MD time step At was found to be dependent on 
the values of f3 and m . In order to simultaneously op- 
timize the acceptance rate, decorrelation and execution 
time, At was adjusted in the range [0.01,0.03], while 
the number of steps N T was chosen randomly from a 
Poisson distribution such that the average MD trajec- 
tory length between updates of the pseudofermion field 
was f = N t At = 2. The choice of f ~ 2.5 was found to 
give optimal decorrelation. 

The HMC algorithm is the method of choice for lat- 
tices larger than 20 2 x 20. As a check on the HMC code, 
the datapoints for 16 2 x 16 computed using the Metropo- 
lis algorithm in Ref. Q were recomputed using the HMC 
method, and found to agree within statistical uncertain- 
ties. In all cases, the uncertainties were estimated using 
the Jackknife method 12911. 



A. The semimetal-insulator transition 

In order to determine the critical coupling f3 c for spon- 
taneous chiral symmetry breaking, we calculated the chi- 
ral condensate a and susceptibility \i f° r P between 0.05 
and 0.5, and for m between 0.0025 and 0.020 (in lattice 
units). Fig. [5] shows our data for lattice sizes 20 2 x 20 
(upper panels) and 28 2 x 28 (lower panels). 

The chiral condensate increases as /3 is decreased, more 
sharply so below (3 ~ 0.1. This behavior becomes more 
pronounced as mo is decreased, providing the first indi- 
cation of a phase transition as the Coulomb coupling is 
increased. In turn, the susceptibility also grows sharply 
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FIG. 3: (Color online) Logarithmic derivative R for lattice 
sizes 20 2 x 20 (left panel) and 28 2 x 28 (right panel). The 
solid lines for the largest four values of /3 correspond to the 
restricted fits shown in Fig. [2] The dashed red line in the right 
panel connects the datapoints for j3 = 1/15.0 s=s 0.067, where 
the downward slope is characteristic of the spontaneously bro- 
ken phase. On the other hand, the dashed blue line connecting 
the datapoints for j3 = 1/11.0 ~ 0.091 clearly indicates that 
chiral symmetry remains unbroken in the limit m — * for 
that value of /3. The evidence for spontaneous chiral symme- 
try breaking is significantly stronger for 28 2 x 28, where the 
data for /3 = 1/13.0 ~ 0.077 are consistent with the broken 
phase, while for 20 2 x 20 the opposite is true. The 28 2 x 28 lat- 
tice favors a slightly larger value of /3 C , while simultaneously 
disfavoring the classical critical exponent 8 — 3. 



FIG. 2: (Color online) Chiral condensate (upper left panel) 
and susceptibility (upper right panel) for lattice size 20 2 x 20. 
Lower panels show the same quantities for 28 2 x 28. The lines 
represent % 2 fits of Eq. (|5.8[) to a only, with X , X 1 ,Y 1 , 8 
and /3 C as free parameters; the datapoints with largest finite- 
size effects have been excluded from the fit. The optimal 
parameter values are: for 20 2 x 20, X = 0.665 ± 0.2, X 1 = 
-0.280 ± 0.088 and Y 1 = -0.2869 ± 0.090, 8 = 2.27 ± 0.13, 
f3 c = 0.0721 ±0.0006; for 28 2 x 28, X = 0.3427 ±0.028, Xi = 
-0.190 ± 0.014 and Y 1 = -0.179 ± 0.014, 8 = 2.309 ± 0.037, 
(3 C — 0.0785 ±0.0003. The uncertainties are purely statistical. 



eludes the use of bare masses m that are small enough 
so that the distortion introduced is negligible. What 
is needed is a controlled way of obtaining information 
about the massless limit, using the data at hand, taken 
at small but finite m . A suitable observable is provided 
by the logarithmic derivative R [3l| of the chiral conden- 
sate with respect to m , 



R 



diner 



d In ; 



da 



(5.1) 



around (3 ~ 0.1. This feature tends to disappear for 
m > 0.010 as the lattice volume is increased. Thus, 
in order to understand the properties of the transition, 
masses smaller than mo ~ 0.010 should be used in the 
simulation. This situation is similar to that encountered 
in quenched QED 4 [3(| where it was concluded that for 
the critical region to be reached, bare masses smaller than 
~ 0.025 should be used. On the other hand, for the small- 
est mass of rn = 0.0025, the change in the susceptibility 
as a function of the lattice volume appears to be rela- 
tively mild for (3 > 0.09. The rise in the susceptibility is 
therefore likely to be a real feature, indicating that the 
critical region has been reached. 

In spite of the compelling qualitative evidence pre- 
sented above, the nature of the simulational study pre- 



which allows for a more precise determination of the crit- 
ical coupling (3 C , as well as for an estimate of the universal 
critical exponent 5 (see Eq. I5.4[) . In the limit m — > 0, 
R — > 1 in the chirally symmetric phase since a oc m ; 
while at the critical coupling = j3 c , one expects R — » 
1 /5. Finally, R vanishes in the limit ra — > in the spon- 
taneously broken phase, where a ^ for m — * 0. The 
data on R in Fig. (right panel) indicate that chiral sym- 
metry is spontaneously broken for (3 = 1/14.0 ss 0.071, 
but remains unbroken for (3 = 1/11.0 » 0.091. We thus 
conclude, using the 28 2 x 28 data, that 

0.071 < (3 C < 0.091, (5.2) 

which could be further refined by use of larger lattice 
volumes and smaller values of m n . 
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B. Determination of the equation of state 

While the logarithmic derivative R may provide model- 
independent information on the critical coupling as well 
as the exponent S, it involves the chiral susceptibility 
and is therefore prone to large finite-size effects. A more 
accurate determination of /3 C can be achieved by means 
of an appropriate equation of state (EOS) 



f(<r,0), 



(5.3) 



which is to be fitted to the simulation data on the chiral 
condensate. This EOS can then yield direct information 
on f3 c as well as the critical exponents 6 and j3, defined 
by 



<91n<7 







d In TO 
d In cr 



0=/3 c ,m o ^o 



In addition, using the scaling relation 

/8(a-i)=7 

one can obtain the critical exponent 7, defined by 
dhix 



7; 



(5.4) 
(5.5) 

(5.6) 
(5.7) 



m =0,/3^/3 c 



The EOS also provides a means for an extrapolation 
?n — * 0, which necessitates an ansatz for Eq. (|5.3p . We 
have considered an EOS similar to those successfully ap- 
plied us, m, m, m to qed 4 , 



(5.8) 



where the functions X and Y are expanded around f3 c 
such that X(fi) = X +X 1 (1 -j3/0 c ) and Y(J3) = Y 1 (l- 
f3/(3 c )- The dependence of Eq. (|5.8p on a is 



A(a) = a\ f 3 (a) = a s 



(5.9) 



where b = 5 — 1 / (5. Thus Eq. (|5.8|) can be used to study 
deviations from the classical exponents (5 = 3 and /3 = 
1/2. It should be noted that for the case of QED 4 [32], 
[33l . [341 . HH , an extended version of the ansatz (|5.9p has 
been used to include logarithmic corrections to the EOS. 

While it is possible to fit both a and \i simultane- 
ously, it is advantageous to use the latter quantity as a 
consistency check only, as the finite-size effects are much 
smaller for a. It is also useful to restrict the fit range 
to the datapoints where such effects are not too large. 
The results of the fits with restricted range are given in 
Fig. [2] whereas the results of a full fit to all datapoints 
is shown in Fig. 2J The results for (3 C and <5 are much 
more consistent for the restricted dataset. The fit re- 
sults for the restricted 28 2 x 28 dataset indicate a critical 
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FIG. 4: (Color online) Chiral condensate (upper left panel) 
and susceptibility (upper right panel) for lattice size 20 x 20. 
Lower panels show the same quantities for 28 x 28. The lines 
represent \ 2 fits of Eq. (|5.8[) to a only, with X ,X 1 , Y 1 ,5 and 
/3 C as free parameters; all datapoints have been included in 
the fit, regardless of the estimated magnitude of finite-size 
effects. The optimal parameter values are: for 20 2 x 20, X = 
0.364 ±0.029, Xi = -0.156 ±0.013 and Fi = -0.159 ±0.013, 
5 = 2.573 ± 0.041, /3 C = 0.0715 ± 0.0003; for 28 2 x 28, X = 
0.834 ±0.024, X 1 = -0.409 ±0.011 and Y 1 = -0.418 ± 0.012, 
S = 1.889 ± 0.017, /3 C = 0.0815 ± 0.0004. 



coupling of j3 c — 0.0785 ± 0.0003 and a critical expo- 
nent 5 = 2.309 ± 0.037. All of the fits described above 
have been performed using the constraint 6=1, which 
is equivalent to the assumption 7=1, using Eq. (|5.6[) . 
However, we have also relaxed this constraint by treating 
b as an additional free parameter in the fit. In all cases, 
no significant deviations from 6=1 were found for any 
of the fits. Nevertheless, it would still be desirable to use 
larger lattices in order to minimize the finite-size effects 
at smaller values of (3. 

However, it is significant that the present results for 
both 20 2 x 20 and 28 2 x 28 favor values of 6 - 2.3 and 
b ~ 1.0, which strongly disfavors the classical mean- field 
exponents 6 = 3, (3 = 1/2. Fits using classical exponents 
tend to become less and less favored when the lattice 
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FIG. 5: (Color online) Fisher plot for lattice sizes of 20 2 x 20 
(left panel) and 28 2 x 28 (right panel). The curved lines con- 
nect datapoints of equal m , from the lower left-hand corner 
these are m = 0.0025, 0.005, 0.010, 0.015 and 0.020. The ver- 
tical lines of equal /3 correspond to the restricted EOS fits 
shown in Fig. [2] The curvature in these lines indicate de- 
viation from the classical critical exponents. At the critical 
coupling, the extrapolation of the lines of equal /3 crosses the 
origin. Finite-size effects tend to turn the lines clockwise. 



volume is increased, which is also reflected in the "Fisher 
plot" shown in Fig. [5] In particular, consistent fits for 5 
can be achieved using data for 20 2 x 20 and 28 2 x 28 if 
the fit range is restricted to those datapoints where the 
finite-size effects are under reasonable control, as shown 
in Fig. H 

It has been argued in Ref. 0] that the semimetal- 
insulator phase transition should present an essential sin- 
gularity, in the sense that the EOS for zero mass in the 
broken phase would be given by 



(7 = C exp ( — 



(5.10) 



with C , C\ constants. This expression has vanishing 
derivatives to all orders at the critical point, and is said 
to be characterized by Miransky scaling [3] ■ The critical 
exponents corresponding to such a transition are 5=1, 
f3 = oo and 7=1. This type of transition has sometimes 
been referred to as a Kosterlitz-Thouless transition, even 
though strictly speaking the latter does not involve spon- 
taneous symmetry breaking. The value 5 = 1 is appar- 
ently ruled out by the considerable dependence of the 
susceptibility on m even for large values of /3, which are 
far from the transition and where the finite-size effects 
are small. If the value of d was close to unity, one would 
observe a susceptibility which is independent of m as 
the critical point is approached. While our data does not 
favor an interpretation in terms of Miransky scaling, a 
full consideration of this issue is beyond the scope of the 
present paper. 



C. Finite-size effects 

If a realistic picture of the properties of the semimetal- 
insulator transition, as exhibited by the low-energy effec- 
tive theory considered here, is to be obtained, a proper 
assessment of the finite-size effects has to be made. In 
general, the lattice volume should ideally be large enough 
such that all explicit degrees of freedom (represented in 
this case by m ) as well as any dynamically generated 
ones (the Goldstone boson associated with spontaneous 
chiral symmetry breaking) can be contained. In order to 
illustrate the finite-size effects, the chiral condensate and 
susceptibility have been plotted for volumes of 20 2 x 20 
and 28 2 x 28 in Fig. 

As expected from the quite different nature of the low- 
energy theory of graphene in the spatial and temporal di- 
rections, the finite-size effects observed in the simulation 
are also different. Increasing the extent of the temporal 
dimension leads to an increase in the condensate er, as 
would be expected by comparison with QED 4 where the 
finite-size effects are dominated by such behavior. The 
finite-size effects in the temporal dimension grow as m 
is decreased, and do not depend strongly on (3. This in- 
dicates that the effects are due to distortion of the stag- 
gered propagator involving the bare mass m . 

On the other hand, increasing the extent L x of the 
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FIG. 6: Chiral condensate (left panel) and susceptibility 
(right panel), shown here for lattices sizes 20 2 x 20 (open 
symbols) and 28 2 x 28 (filled symbols). The finite -size effects 
are typically much larger for the susceptibility, and become 
large for the condensate as well at small /3, particularly in the 
broken phase. Errors shown are purely statistical, and are in 
most cases comparable to the size of the datapoints. 
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spatial directions has a quite different effect on the chiral 
condensate. The effect is to lower the value of a, which 
is opposite to the effect of increasing L t . The relative 
change in a also appears to be roughly independent of 
m , such that the absolute shift is larger for larger values 
of m . It is also noteworthy that the finite-size effects 
in L x are very small in the unbroken phase, as shown in 
Fig. while they quickly become large with the onset 
of spontaneous chiral symmetry breaking. We therefore 
conclude that these effects are due to the emergence of a 
dynamically generated Goldstone mode which is highly 
spatially extended. 



For all the results presented here, the extent L z — 8 has 
been used for the bulk dimension, in which the fcrmionic 
degrees of freedom do not propagate. Increasing the size 
of that dimension has apparently no effect whatsoever 
on the results for the chiral condensate and susceptibil- 
ity, as demonstrated by a comparison between results on 
a 14 2 x 14 lattice with L z — 8 and L z = 14. The results 
for all observables in question are compatible within sta- 
tistical uncertainties, and binning of the data for a into a 
histogram plot also shows no perceptible differences be- 
tween the two event distributions. Apparently, restrict- 
ing the bulk dimension to L z — 8 has no significant effect 
on the accuracy of our results, although an increased L z 
can be accommodated if necessary, as this has little ef- 
fect on the total computational cost. Such a result is 
nevertheless somewhat intuitive, as the fermions do not 
propagate in the bulk, and thus should be mostly insen- 
sitive to the presence of a boundary in that dimension. 
However, it should still be noted that [33, [33|, [34j in the 
context of QED 4 the main effect of the boundary is to in- 
troduce a constant background component into the gauge 
fields. In other words, at finite volume the results can be 
well described in terms of a renormalized staggered lat- 
tice propagator, augmented by a constant background 
field that may vary from one configuration to the next. 



In addition to shifting the calculated values of the con- 
densate, finite-size effects may also influence the distri- 
bution of the measured MC samples. We have observed 
that for small lattice volumes, the simulation exhibits a 
tendency to jump between two different states, akin to 
the effect noted in the QCD simulations of Ref. 24]. This 
effect appears to be strongest in the quenched case, and 
weakens as more fermion flavors are added. The area of 
parameter space most affected is just above C , where 
the Coulomb interaction is not yet quite strong enough 
to break the chiral symmetry, and a is strongly fluctu- 
ating. As this first-order feature also tends to disappear 
with increasing decorrelation and decreasing finite-size 
effects, we attribute it to a combination of these factors. 
This is in line with Ref. (24[, where attempts to fit the 
event distribution with two Gaussians did not turn out 
satisfactorily. 



VI. TESTS AND CROSS-CHECKS 

In this section, we briefly describe the various tests 
performed in order to validate our simulations. Using the 
formalism described in Sec. IIII[ we extended our code to 
perform simulations of QED in 2+1 dimensions (QED 3 ), 
and compared our results with those from Ref. [30(. In 
this case the differences with graphene are that the gauge 
field lives in one less spatial dimension, and that all 
the components of the gauge field are dynamical, since 
Lorentz invariance is respected. 

We have also developed another test based on QED in 
3+1 dimensions (QED/), which we compared with the 
results of Ref. [32|, |33, ]34 1 . In this case the differences 
with graphene affect the fermion field, which lives in one 
more dimension. As in the previous case all the com- 
ponents of the gauge field contribute, as the theory is 
Lorentz invariant. Our lattice Monte Carlo implementa- 
tion has satisfactorily passed all of the abovementioncd 



QED 3 , 8 d (Kogut) N f = 
8 3 (Kogut) N f : 
8 3 (Kogut) N f : 
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FIG. 7: Data of Ref. [301 ] compared with our implementation 
of the QED 3 simulation. The filled datapoints are our results, 
whereas the empty ones denote the results of Ref. [3uT ]. The 
lines connecting the datapoints of Ref. [3(j are intended as a 
guide to the eye. 
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tests. A comparison between our results for QED 3 and 
those of Ref. [30] is shown in Fig. [7J 

In addition to these major checks, the following usually 
overlooked ones were also performed: explicit verification 
using a Computer Algebra System (CAS) of the correct 
structure of the staggered fermion operator, invariancc 
of the action and the observables under gauge transfor- 
mations, and reversibility of the HMC algorithm within 
each MD evolution. We finally note that the computing 
time required by the present calculations is ~ 10 5 CPU- 
hours, which is in line with an estimate given by Hands 
and Strouthos in Ref. Allocations of this size are 

routinely available at various supcrcomputing centers. 



VII. OBSERVATION OF THE TRANSITION 

The experimental detection of excitonic instabilities 
in graphene depends on the size of the induced gap A. 
Unfortunately, computing A in absolute units requires 
knowledge of a suitable dimensionful observable (other 
than A itself) to calibrate the calculation. To our knowl- 
edge such a quantity is not yet available. In Ref. 0], the 
excitonic gap was estimated within a gap-equation ap- 
proach by assuming a value of the cutoff of the order of 
the inverse lattice constant of the graphene honeycomb 
lattice. In that study, the gap was found to be of the 
order of a few tens of K. However, such a procedure only 
constitutes an order-of-magnitude estimate. As in our 
approach, the size of the gap can be determined in abso- 
lute units only after calibration of the calculation using 
a dimensionful observable. 

Another issue of significance from the experimental 
point of view is the effect of impurities and lattice de- 
fects. These were investigated in Ref. (36|, and they were 
found to have a substantial impact on the low-energy ex- 
citations in graphene. Also, Ref. [13] has studied the sta- 
bility of the excitonic insulating phase in the presence of 
impurities, lattice defects and thermal fluctuations, and 
concluded that all of these effects tend to suppress the 
excitonic instability. Clearly, the experimental demon- 
stration of the semimetal-insulator transition in graphene 
will be challenging from the point of view of sample qual- 
ity. 

As the mere presence of a substrate will likely eliminate 
the insulating phase due to screening of the Coulomb in- 
teraction, the most favorable experimental setup would 
involve samples of suspended graphene. Fortunately, this 
may also serve to eliminate most of the abovementioned 
concerns. Indeed, it has recently been found in Ref. [3] 
that in order to access the intrinsic electronic properties 
of graphene, thorough current-annealing of suspended 
samples is necessary. The annealed samples were found 
to exhibit a greatly improved carrier mobility, far in 
excess of the values reported for conventional samples 
on a substrate. Also, the demonstration of Shubnikov- 
deHaas (SdH) oscillations suggests that the mean free 
path in current state-of-the-art suspended graphene is 



comparable to presently achievable sample dimensions of 
a few /im. Thus, graphene samples of sufficient qual- 
ity to demonstrate the excitonic instability will likely be 
available in the near future. 

To summarize, our work in Ref. [5[ indicates that the 
excitonic insulating effect in graphene is unlikely to be 
observed unless the graphene sheet is freely suspended, 
such that the Coulomb interaction is not screened by the 
dielectric substrate. Further, the experimental work in 
Ref. 0] has demonstrated that the elimination of impuri- 
ties and defects is necessary in order to access the intrin- 
sic electronic properties of graphene. As both of these 
conditions can nowadays be fulfilled by experiment, we 
hope that the appearance of the excitonic gap will be 
demonstrated in the near future. 



VIII. CONCLUSIONS 

We have described the low-energy effective theory of 
graphene, its gauge and global symmetries, and shown 
how a discretized lattice formulation can be constructed 
such that it contains the correct number of degrees of 
freedom and partially retains chiral invariance at finite 
lattice spacing. We have also explained in detail the 
numerical methods employed to perform lattice Monte 
Carlo simulations of the discretized theory, focusing on 
the determination of the location and properties of the 
semimetal-insulator phase transition. 

On the theoretical side, we conclude that our extended 
analysis is consistent with the findings of Ref. [f| , which 
predict that suspended graphene should possess an ex- 
citonic gap in the band structure. We have now, using 
the HMC algorithm, extended the results of Ref. [5[ to 
much larger lattice volumes, as well as smaller fermion 
masses. While the scenario first reported in Ref. [5[ is 
confirmed by the present results, the larger lattices used 
also provide tantalizing hints that the phase transition 
is not of infinite order, as predicted in Ref. 0, nor is it 
likely to be described by classical critical exponents. In 
order to achieve a precise determination of the critical 
exponents it is necessary to perform simulations at much 
larger lattices, potentially as large as 48 2 x 48. We are 
currently exploring the feasibility of such simulations by 
benchmarking our code on a 36 2 x 36 lattice. 

An accurate determination of the critical coupling and 
the critical exponents will provide a solid understanding 
of the universality class of this transition, as well as an- 
other piece of experimentally verifiable information on 
the electronic properties of graphene. 
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